Curso I – SISSA

Estadísticas Espaciales y Temporales

Objetivo y contenido

Objetivo y contenido

Objetivo: realizar estadísticas espacio-temporales básicas utilizando datos grillados y vectoriales.

  • Estadísticas ráster
    • Mínimo, máximo, media, desviación estándar, y suma
    • Frecuencia de celdas en un ráster
    • Exportar valores derivados de un ráster como una serie temporal
    • Ejemplo: estadísticas ráster sobre un área
    • Calculando la media de P y Ea sobre un área
    • Analizando patrones de P - Ea

Estadísticas ráster

Estadísticas ráster

Comencemos con algunas estadísticas ráster básicas. Para esto, necesitamos cargar el paquete terra en R:

library(terra)

Por favor, cargue el siguiente ráster de temperaturas promedio mensuales globales para el mes de febrero.

temp_feb_path <- "C:/Users/Temperature_Data/wc2.1_2.5m_tavg_02.tif"
temp_feb      <- rast(temp_feb_path)

Nota: los datos se proporcionan en grados Celsius.

Estadísticas ráster

Empecemos por mirar los datos:

plot(temp_feb)

Estadísticas ráster

hist(temp_feb)

Valores mínimos y máximos

Podemos extraer los valores mínimos y máximos de un ráster:

minmax(temp_feb)
##     wc2.1_2.5m_tavg_02
## min             -44.80
## max              32.96

Valores mínimos y máximos

Si solo queremos verificar los valores mínimos y máximos (es decir, no guardar el número ni usarlo en un cálculo), los valores aparecen cuando escribimos el nombre del objeto y presionamos enter.

temp_feb
## class       : SpatRaster 
## dimensions  : 4320, 8640, 1  (nrow, ncol, nlyr)
## resolution  : 0.04166667, 0.04166667  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : wc2.1_2.5m_tavg_02.tif 
## name        : wc2.1_2.5m_tavg_02 
## min value   :             -44.80 
## max value   :              32.96

Celdas con atributos específicos

Para averiguar el número de celdas que tienen un cierto valor, podemos usar la función freq.

freq(x, digits = 0, value = NULL, ...)

x: capa ráster

digits: se utiliza para redondear los valores de las celdas

value: valor opcional. Cuenta solo las celdas con este valor

Celdas con atributos específicos

Si no especificamos un valor o una condición de redondeo:

freq(temp_feb)
##    layer value  count
## 1      1   -45  17478
## 2      1   -44  83878
## 3      1   -43 117982
## 4      1   -42 102845
## 5      1   -41 106076
## 6      1   -40 152928
## 7      1   -39 154018
## 8      1   -38 151494
## 9      1   -37 226850
## 10     1   -36 213282
## 11     1   -35 232837
## 12     1   -34 293376
## 13     1   -33 273896
## 14     1   -32 275022
## 15     1   -31 241752
## 16     1   -30 213349
## 17     1   -29 250686
## 18     1   -28 211607
## 19     1   -27 203796
## 20     1   -26 198331
## 21     1   -25 219682
## 22     1   -24 205186
## 23     1   -23 215089
## 24     1   -22 186505
## 25     1   -21 182648
## 26     1   -20 197982
## 27     1   -19 194459
## 28     1   -18 210066
## 29     1   -17 241606
## 30     1   -16 259705
## 31     1   -15 223245
## 32     1   -14 209189
## 33     1   -13 211503
## 34     1   -12 191229
## 35     1   -11 157554
## 36     1   -10 149361
## 37     1    -9 136770
## 38     1    -8 130283
## 39     1    -7 125210
## 40     1    -6 121743
## 41     1    -5 110963
## 42     1    -4 109459
## 43     1    -3 116599
## 44     1    -2 111883
## 45     1    -1  94808
## 46     1     0 117162
## 47     1     1  94409
## 48     1     2  84865
## 49     1     3  77204
## 50     1     4  84120
## 51     1     5  83697
## 52     1     6  80351
## 53     1     7  74984
## 54     1     8  72369
## 55     1     9  73109
## 56     1    10  73799
## 57     1    11  67978
## 58     1    12  84051
## 59     1    13 101461
## 60     1    14 119884
## 61     1    15 110814
## 62     1    16  88892
## 63     1    17  85445
## 64     1    18 104649
## 65     1    19 121411
## 66     1    20 140826
## 67     1    21 172560
## 68     1    22 214183
## 69     1    23 237690
## 70     1    24 287524
## 71     1    25 399249
## 72     1    26 463549
## 73     1    27 335760
## 74     1    28 192903
## 75     1    29 104746
## 76     1    30  82156
## 77     1    31  50880
## 78     1    32  15303
## 79     1    33    419

Celdas con atributos específicos

Ahora, vamos a buscar solo las celdas con NA:

freq(temp_feb, value = NA)
##   layer value    count
## 1     1    NA 24792188

Ahora, busquemos el valor 20 grados, pero dentro de un margen de redondeo de 0.1 (ten en cuenta que esto significa que solo se incluyen los valores de 19.95 a 20.05).

freq(temp_feb, digits = 1, value = 20)
##   layer value count
## 1     1    20 13606

Exportar series temporales derivadas de un ráster

Utilicemos lo que acabamos de aprender para exportar una serie temporal de la temperatura media mensual global durante los 12 meses para los que tenemos datos (sobre áreas terrestres, ya que los océanos se almacenan como NA).

Comencemos cargando los rásters y creando un stack con las 12 capas:

t_path <- "C:/Users/Temperature"
t_files <- list.files(t_path, full.names = TRUE)
t_stack <- rast(t_files)

Exportar series temporales derivadas de un ráster

Ahora queremos calcular el valor medio en cada capa:

mean_values <- global(t_stack, stat = mean, na.rm = TRUE)
print(mean_values)
##                          mean
## wc2.1_2.5m_tavg_01 -6.9901310
## wc2.1_2.5m_tavg_02 -8.2640013
## wc2.1_2.5m_tavg_03 -8.4594775
## wc2.1_2.5m_tavg_04 -6.3245406
## wc2.1_2.5m_tavg_05 -3.3703036
## wc2.1_2.5m_tavg_06 -0.6429807
## wc2.1_2.5m_tavg_07  0.1324651
## wc2.1_2.5m_tavg_08 -0.6111403
## wc2.1_2.5m_tavg_09 -2.2624299
## wc2.1_2.5m_tavg_10 -4.0777166
## wc2.1_2.5m_tavg_11 -5.3751544
## wc2.1_2.5m_tavg_12 -5.8202356

Exportar series temporales derivadas de un ráster

Pregunta: Debido al sistema de proyección de los datos en este ejemplo (WGS84), ¿los resultados serán representativos de las medias de temperatura globales (sobre la tierra) o estarán sesgados? Si están sesgados, ¿de qué manera (es decir, son los resultados de nuestro análisis más altos o más bajos que la temperatura promedio real sobre la tierra)?

Exportar series temporales derivadas de un ráster

Veamos los resultados.

plot(1:12, mean_values$mean, type = "l", xaxt = "n", xlab = "Month", ylab = "Temp [deg. C]")
axis(1, 1:12, month.abb)

xaxt = “n” es para excluir las etiquetas en el eje x. Esto se debe a que queremos agregar las abreviaturas de los meses (vea la siguiente línea) en lugar de un valor predeterminado de 1 a 12.

month.abb es un objeto predefinido en R que tiene las abreviaturas de cada mes.

Exportar series temporales derivadas de un ráster

Exportar series temporales derivadas de un ráster

Finalmente, queremos exportar estos valores como una serie temporal.

write.table(mean_values, "C:/User/Exported_TS.csv", row.names = FALSE, sep = ",")

Resumen del proceso

Acabamos de mostrar cómo calcular el valor medio (u otra variable) para cada paso de tiempo. Luego, generamos una serie temporal para los 12 meses – con un valor para cada mes –.

Próximo paso

Ahora queremos calcular la temperatura media anual en cada celda de la grilla. Calcularemos para cada celda la media de los 12 valores que corresponden a los 12 meses.

Calculando valores medios en cada celda

Utilizaremos los mismos datos de temperatura que utilizamos en el ejemplo anterior. Esta vez, en lugar de global, podemos utilizar la función mean.

mean_raster <- mean(t_stack)
plot(mean_raster)

Calculando valores medios en cada celda

Visualización de valores NA

A veces es útil visualizar la ubicación de los valores NA en un archivo ráster. Para este propósito, podemos usar el parámetro colNA de la función plot.

plot(mean_raster, colNA = "red")

Ejemplo: Estadísticas ráster (P - Ea) sobre un área

Estadísticas ráster

Para este ejercicio, queremos analizar los patrones de precipitación y evapotranspiración en Uganda. Tenemos los siguientes datos:

  • Shapefile para los países de la cuenca del Nilo
  • P mensual global de CHIRPS para 2018
  • Ea mensual sobre África de SSEBop para 2018

Queremos generar un gráfico de línea para la media mensual de P y la media mensual de Ea en Uganda.

Estadísticas ráster

Paso 1: Cargar los paquetes necesarios en R.

library(sf)
library(terra)

Paso 2: Establecer las rutas para el shapefile y las carpetas con los datos de CHIRPS y SSEBop.

dir <- "C:/Users/L5_Spatial_and_Temporal_Statistics_Data"
shape_path  <- file.path(dir, "Shapefile/Nile_Basin_Countries_GAUL2014_2.shp")
chirps_path <- file.path(dir, "CHIRPS")
ssebop_path <- file.path(dir, "SSEBop")

Estadísticas ráster

Paso 3: Leer el shapefile y extraer el polígono correspondiente a Uganda.

countries <- read_sf(shape_path)
uganda    <- countries[5,]
uganda <- vect(uganda)

Estadísticas ráster

plot(uganda)

Estadísticas ráster

Paso 4: Cree un stack de los archivos mensuales de CHIRPS.

chirps_files <- list.files(chirps_path, full.names = TRUE)
chirps       <- rast(chirps_files)
plot(chirps)

Estadísticas ráster

Queremos calcular la media mensual de precipitación sobre Uganda. Eso significa que solo queremos calcular estadísticas ráster dentro de este shapefile. ¿Cuál la forma de hacer esto?

  • Enmáscar el ráster con el shapefile (¡no recortar!)

  • Calcular las estadísticas utilizando este ráster enmascarado

Pero, si buscamos otras opciones, podemos ver que existe una función que ya realiza este proceso por nosotros: la función extract (en el paquete terra).

Estadísticas ráster

Paso 5: Extraiga el promedio mensual de P sobre Uganda.

ug_pmean <- terra::extract(chirps, uganda, fun = mean, na.rm = TRUE)

Podemos ver que nuestra salida está almacenada como una matriz. El primer valor debe ser eliminado.

print(ug_pmean)
##   ID chirps_monthly2018.01 chirps_monthly2018.02 chirps_monthly2018.03
## 1  1              21.60957              46.26966               171.186
##   chirps_monthly2018.04 chirps_monthly2018.05 chirps_monthly2018.06
## 1              227.9121              187.1596              125.8768
##   chirps_monthly2018.07 chirps_monthly2018.08 chirps_monthly2018.09
## 1              62.84645              134.0909              103.6331
##   chirps_monthly2018.10 chirps_monthly2018.11 chirps_monthly2018.12
## 1              126.2652              86.03881              99.42722

Estadísticas ráster

Convirtamos la salida en un vector numérico para que sea más fácil de manejar.

ug_pmean <- as.numeric(ug_pmean)[-1]
print(ug_pmean)
##  [1]  21.60957  46.26966 171.18602 227.91208 187.15955 125.87678  62.84645
##  [8] 134.09085 103.63311 126.26518  86.03881  99.42722

Estadísticas ráster

Paso 6: Repita el mismo proceso para los datos de Ea: primero hay que crear un ráster stack y después extraer el promedio mensual de Ea.

ssebop_files <- list.files(ssebop_path, full.names = TRUE)
ssebop       <- rast(ssebop_files)
plot(ssebop)

Estadísticas ráster

Estadísticas ráster

ug_eamean <- terra::extract(ssebop, uganda, fun = mean, na.rm = TRUE)
ug_eamean <- as.numeric(ug_eamean)[-1]
print(ug_eamean)
##  [1]  71.16544  50.92773  93.04707  95.88938  98.27621  93.29266  85.61035
##  [8]  87.66707  96.88887 100.74566  94.00632  91.06700

Estadísticas ráster

Paso 7: Cree un gráfico de líneas para el promedio mensual de P y Ea sobre Uganda.

Primero, comenzamos graficando solo el promedio mensual de P:

plot(1:12, ug_pmean, type = "l", main = "Promedio mensual de P y Ea para Uganda - 2018",
col = "blue", lwd = 2, xlab = "Mes", ylab = "P y Ea [mm]", xaxt = "n")

Luego, podemos agregar el Ea usando el comando lines.

lines(1:12, ug_eamean, col = "darkred", lwd = 2)

Estadísticas ráster

Finalmente, podemos incluir las abreviaturas de los meses en el eje x, las líneas de la cuadrícula y una leyenda.

axis(1, 1:12, month.abb)
grid()
legend("topright", c("P", "Ea"), lty = 1, col = c("blue", "darkred"))

Estadísticas ráster

Estadísticas ráster

Para este ejercicio, generaremos tres salidas más:

  1. Boxplot de la evaporación mensual en Uganda

  2. Distribución espacial de la precipitación anual acumulada menos la evaporación en Uganda

  3. Histograma de los valores de la precipitación anual acumulada menos evaporación en Uganda

Estadísticas ráster

Para poder completar nuestros tres objetivos adicionales, queremos crear rásters de P y Ea mensuales enmascarados para Uganda.

Use la función mask para generar los rásters de P y Ea mensuales.

uganda_prec <- crop(chirps, uganda, snap = "out")
uganda_prec <- mask(uganda_prec, uganda)
uganda_ea <- crop(ssebop, uganda, snap = "out")
uganda_ea <- mask(uganda_ea, uganda)

Estadísticas ráster

plot(uganda_prec)

Estadísticas ráster

plot(uganda_ea)

Estadísticas ráster - a)

Generar un boxplot para mostrar la dispersión de los valores de Ea para cada mes.

boxplot(uganda_ea, main = "Evaporación mensual en Uganda - 2018", col = "red", xlab = "Mes", ylab = "Ea [mm]", xaxt = "n")
axis(1, 1:12, month.abb)
grid()

Estadísticas ráster - b)

Para calcular P menos Ea para 2018, primero debemos calcular la precipitación y la evaporación anuales acumuladas para Uganda.

Generar rásters para la precipitación y la evaporación anuales acumuladas:

annual_prec <- sum(uganda_prec)
annual_ea   <- sum(uganda_ea)

Estadísticas ráster - b)

plot(annual_prec, main = "P anual para 2018")

Estadísticas ráster - b)

plot(annual_ea, main = "Ea anual para 2018")

Estadísticas ráster - b)

Remuestrear el ráster de P a la misma resolución espacial que el ráster de Ea.

annual_prec <- resample(annual_prec, annual_ea, method = "bilinear")

Calcular la precipitación anual acumulada menos la evaporación en Uganda.

p_minus_ea <- annual_prec - annual_ea

Estadísticas ráster - b)

plot(p_minus_ea, main = "P - Ea para 2018")

Estadísticas ráster - c)

Realizar un histograma de la distribución de los valores de la precipitación anual acumulada menos la evaporación.

hist(p_minus_ea)

Graficando Shapefiles y rásters

Recuerda, si te sientes mejor graficando shapefiles y rásters utilizando otro programa (por ejemplo, ArgGIS, QGIS, Google Earth Engine), esto aún puede hacerse fácilmente. Solo recuerda que primero debes guardar tus resultados (es decir, exportar los datos).

¡Gracias por su atención!